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ABSTRACT 

We report new methods for evaluating realistic observing programs that 
search stars for planets by direct imaging, where observations are selected from 
an optimized star list, and where stars can be observed multiple times. We show 
how these methods bring critical insight into the design of the mission and its 
instruments. These methods provide an estimate of the outcome of the observing 
program: the probability distribution of discoveries (detection and/or character- 
ization), and an estimate of the occurrence rate of planets (77). We show that 
these parameters can be accurately estimated from a single mission simulation, 
without the need for a complete Monte Carlo mission simulation, and we prove 
the accuracy of this new approach. Our methods provide the tools to define a 
mission for a particular science goal, for example defined by the expected num- 
ber of discoveries and its confidence level. We detail how an optimized star list 
can be built and how successive observations can be selected. Our approach also 
provides other critical mission attributes, such as the number of stars expected 
to be searched, and the probabihty of zero discoveries. Because these attributes 
depend strongly on the mission scale (telescope diameter, observing capabilities 
and constraints, mission lifetime, etc.), our methods are directly applicable to 
the design of such future missions and provide guidance to the mission and in- 
strument design based on scientific performance. We illustrate our new methods 
with practical calculations and exploratory design reference missions (DRMs) for 
the James Webb Space Telescope (JWST) operating with a distant starshade 
to reduce scattered and diffracted starlight on the focal plane. We estimate 
that 5 habitable Earth-mass planets would be discovered and characterized with 
spectroscopy, with a probability of zero discoveries of 0.004, assuming a small 
fraction of JWST observing time (7%), rj = 0.3, and 70 observing visits, limited 
by starshade fuel. 
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Introduction 



Various proposals for direct-imaging searches for Earth-like planets are now on the table 
dCuvon et aPboOsi ICash et al.lboosi : ISpergel et aPboOol : ISoummer et al.lboOQbh . This paper 
presents methods to guide the development of the new instrumentation that such searches 
will require. They can also help optimize observing programs by scheduling observations to 
maximize their impact. Using these new methods, we can predict the outcome of a single 
observation or an entire observing program, and interpret the observational results. These 
tools should be widely useful for studying, comparing, and optimizing alternative direct- 
search concepts. 

The completeness of a direct limiting search observation (LSO) is the fraction of all 
possible planets of interest (POIs) that satisfy the detection criteria. By definition, an 
LSO has sufficient exposure time, texpj' on the j*^ star, to reach the systematic limit of 
the instrument. For optical missions, an LSO is assumed to achieve the desired photometric 
signal-to-noise ratio on the limiting source, which has magnitude magj + AmagQ, where magj 
is the stellar magnitude, and AmagQ is the limiting magnitu de difference with the star (flux 
contrast), expected to be determined by speckle instability (IBrownl 120051 ). 



In the simplest picture, an LSO discovers any and all POIs that satisfy two criteria at 
the time of the observation: 

Amag < AmagQ , (1) 

and 

s > IWA , (2) 

where s is the angular separation between planet and star, and IWA is the inner work- 
ing angle (angular radius of the real or effective central field obscuration). This simplified 
picture of a sharp dividing line between detectable and undetectable planets has proven 
useful for measuring search power to first order. Fidelity could be improved, if necessary, 
by adding more complex detection criteria, such as a detection probability that varies in a 
more complicated and realistic way over the field of view. 

In a common treatment, which we follow here, the POIs are body-twins of Earth with 
orbits in the habitable zone. 

Let a fraction t] of all stars in the universe have a POL The list of candidate stars for 
the observing program is a subset of all stars, and within this subset, the expected fraction 
is rj, but the actual fraction will vary. Also, we have limited knowledge of 77, a number 
that is both interesting scientifically and useful operationally to schedule observations to 
maximize discoveries. Therefore, we will draw distinctions when necessary between the true 
value (r^true), which is unknown, the assumed value for science operations (^^ops), which is 
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used to estimate discovery probabilities, and the value estimated from the results of a search 
program [E{ri)]. 

Our new rnethods ext end the original concept of direct-search completeness for exo- 
planets ( iBrownl |2004J . 120051 ) for purposes of optimizing the timing of revisits to previously 
searched stars, increasing the realism of instrument compari sons, and provi ding an estimate 
of 7] from the results of observing programs. Meanwhile, iBrownl (l2009al Jbl) has extended 
completeness studies to indirect detection by reflex astrometry, photometric detection in the 
case of no occultations, and the estimation of orbital parameters from Keplerian data sets. 



2. Evolution and Estimation of Completeness 

2.1. Four Types of Completeness 

Four types of completeness pertain to a program of LSOs: virgin, dynamic, accumulated, 
and ultimate. 

Until now, the scheduling of LSOs in mission studies of direct-search power has been 
based on "virgin," time-i ndependent, first- visit completeness Ci^,-, which is the comp leteness 
of the first ISO of star j JSrownlbool : kgollboO?! : ISavranskv. Kasdin. fc Cadvlboiol ). In this 
study, we include the possibility of multiple LSOs of any target star. 

A non-detection by an LSO rules out some fraction of possible planets: those with 
sufficient angular separation and brightness at the time t of the LSO to be detectable — if 
they existed. After the LSO, a pool of possible planets may remain, comprising planets 
that had not been ruled out by previous LSOs and also were not detectable by the most 
recent LSO. As time goes on, each planet in this pool moves along its unique orbit and may 
become detectable at some future time. In this way, the fraction of all possible planets that 
is detectable on the i*^ visit to the j^^ star at time t — which is the dynamic completeness, 
Cij{t) — depends on the elapsed time since each of the i — 1 preceding LSOs. 

Accumulated completeness (Cjj) is the sum-total completeness of i LSOs of star j: 



a 



(3) 



1=1 



Cij increases monotonically with i. 

Ultimate completeness (Cooj) is the maximum value of Cij. It is the value of complete- 
ness that would be accumulated from an arbitrarily high number of LSOs spread arbitrarily 
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over time: 

oo 

Coo,, = 5^ Q,, (to . (4) 

1=1 

CoD,j < 1 whenever some POI orbits are permanently fainter than AmagQ or permanently 
obscured inside IWA — or never brighter than AmagQ and located outside IWA at the same 
time. 



2.2. Estimating Dynamic Completeness Cj 

All types of completeness are derived from q^-, and all results (probability of discovery 
and mission outcomes) are ultimately based on this quantity. 

Because the planetary position is determined by a transcendental equation (Kepler's 
Equation), must be estimated by Monte Carlo trials. In these trials, we represent the 
universe of POIs by a large random sample of A^^o particular POIs, each of which is defined 
by randomly chosen values for ten parameters: {a, e, Mq, i, Wp, fi, T, -Rp, g, $}, where a 
is the semimajor axis, e is the orbital eccentricity, Mq is the mean anomaly at some definite 
time, i is the inclination angle, Wp is the argument of periastron, f2 is the position angle of 
the ascending node, T is the orbital period, i?p is the planetary radius, q is the geometric 
albedo, and $ comprises the necessary sub-parameters for defining the phase function. 

In this paper, the POIs are Earth-twins on habitable-zone orbits. The particular values 
of six parameters are drawn from random deviates: 0.7a/L < a < LSv^L, < e < 0.35, and 
< Mo < 277 (uniformly distributed; L is the stellar luminosity); i, Wp, and VL uniformly 
distributed on the sphere. Three parameters are delta functions: -Rp = i?©, q depends on 
the filter passband, and $ is the Lambertian phase function. T is given by a and the stellar 
mass via Kepler's Third Law. 

When estimating cij-, for the first LSO, Amag and s are computed for all Nq POIs at, 
say, if: = 0, when the orbital phases of all POIs in the sample are equal to {Mq}. Thereafter, 
record is kept of the epoch of each LSO, ti>i, only for the POIs that have not yet been 
eliminated and are still in play. 

For the i*^ LSO, we identify and count the number of POIs that satisfy Eqs. (1-2), Nij, 
producing: 



Based on Eq. ([5]), all types of completeness can be computed with any required precision by 
appropriately choosing the value of Nq. 
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2.3. Rebound of Cij 



Ci+ij{t) rebounds following the z*^ LSO. Ci+ij{ti+e) = 0, where e is a diminishingly small 
increment of time, and then it rebounds towards a constant, asymptotic value, Ci+ij{oo), as 
the orbits of still-possible POIs lose orbital phase coherence. 

The blue points in Figure 1 show how C2,t rebounds after a first LSO a. t t = of 



HIP 29271 for IWA = 0.075 arcsec ( TH E! A: \S&vranskv. Kasdin. fc Cadvll2010f l q = 0.26 



and AmagQ = 26 (typical value; iBrownl l2005l ). The rebound resembles the response of 



the output voltage for an underdamped, series, LRC circuit after a step change in the 
input voltage. After a linear rise (curved on a logarithmic plot), C2j undergoes damped 
oscillations. The details — rise time, damping time, and asymptotic value C2j(oo) — depend 
on the star's particular mix of habitable-zone orbits resolved by IWA. For different stellar 
mass, luminosity, and distance, for other values of IWA and AmagQ, and for other definitions 
of POIs, the rebound of C2j will be qualitatively similar to Figure 1, but with different details 
and numbers. For example, other factors being equal, a star at larger distance would have a 
longer rebound time, because the central field obscuration (IWA) would limit the observation 
to planets with wider separations, and therefore a longer time would be necessary to loose 
orbit coherence. 

When the next observation is being planned, we need to know Cij{t) for each star in 
play, and where i — 1 observations of star j have already been performed. When i = 1, 
we would use tabulated values of time-independent virgin completeness, and no real-time 
computation would be required. When z > 1, however, we may need an efficient function or 
procedure for computing the values of Cij{t). 

The most accurate, brute-force method would perform a blue-point-type calculation (see 
Figure 1) for every star in play every time a new observation is planned. The number of 
times would be of order the number of stars times the number of observations. For example, 
the number of blue-point-type calculations would exceed 10^ for a program of 100 stars and 
1,000 LSOs, typical for a 4-m class instrument with IWA = 0.075 arcsec. Monte Carlo full- 
mission studies would be impractical, as each of the 400 blue points in Figure 1 took ~5 sec 
to compute on a 3 GHz Intel Xenon processor running Mathematica 6. Therefore, we 
must look at two approximate functions for Cjj(t), one of which may be perfectly adequate 
for first-order scheduling studies. They demand only one or four blue-point-type calculations 
performed a number of times that is of order the number of observations. 

The first alternative approximation is a linear function, Cij{t) = slope x t for t < 
breaktime, and Cij{t) = Cjj(oo) for t > breaktime. This is illustrated by the red curve in 
Figure 1. 
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Fig. 1. — Typical rebound of dynamic completeness C2,j{t) following the first limiting search 
observation (LSO) at t = 0. Blue: values calculated from Eq. 5 with Nq = 20, 000, for Earth- 
like, habitable-zone planets around HIP 29271, assuming IWA = 0.075 arcsec, q = 0.26, and 
AmagQ = 26. Red: the linear approximation. Here, C2j(oo) = 0.18 and breaktime = 
5 X 10'^ sec. 
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We used the following algorithm to find the parameters slope and breaktime = Ci j{oo) / slope. 
The algorithm comprises four blue-point-type computations of Ci .j(t). First, we estimate 
Qj(oo) by computing Cjj(10"'^'' sec). Second, we compute Cjj(10^'^ sec), and use it to make a 
first estimate of the slope, slopel = Ci j{10^-^ sec)/10^'^ sec. (The starting point t = 10^-^ is 
somewhat arbitrary. It should be large enough to afford an accurate value of according 
the counting statistics of Eq. (5), but also confidently smaller than the true value of break- 
time.) Third, we compute a first estimate of breaktime, breaktimel — Cij(10^° sec)/s/o^»el, 
and compute Cij {breaktimel), which produce second estimates of the breaktime and slope: 

Ci A breaktimel) — Cii (10^ sec) 

breaktime2 = breaktimel - ^ , (6) 

slopel 

and slope2 = Cij{10^^) / breaktime2. Fourth, we compute Cij{breaktime2) — the fourth and 
last blue-point-type computation — and use it to produce the final estimates 

^ , , , ^ Ci j(breaktime2) — Ci j( breaktimel) 
breaktimeFmal = breaktime2 — ; — , (7) 

slope2 ' ^ ' 

and slopeFinal = Cj ,,(10^° sec) / breaktimeFinal. (The linear function required about ~20 sec 
to calculate on a 3 GHz Intel Xenon processor running Mathmatica 6.) 

The second alternative approximation is a step function: Cij{t) = for t < breaktime, 
and Cijit) = Ci j{oo) for t > breaktime, where breaktime ~ 10^ sec for habitable-zone orbits 
and popular instrument concepts. For first-order investigations, accurate knowledge of the 
completeness rebound may not be important to the outcome. The important thing is to avoid 
the mistake that breaktime is zero or too small, which error may cause bogus observations 
to pile up on the high-completeness, low-exposure-time stars. (The step function required 
~5 sec to calculate on a 3 GHz Intel Xenon processor running Mathmatica 6.) 

In Figure 2, the first (linear) approximation of Cij{t) is used to illustrate a hypothetical 
program of five LSOs to HIP 29271 starting at absolute time t — 10^ sec. (The abscissa 
here is now linear.) Such approximations would be used only in the scheduling process of 
a DRM, where it may be necessary to compute dynamic completeness for many stars on 
the fly. After the decision observe a particular star, Eq. (5) would be used to produce an 
accurate value for the record. Note that as the number of LSOs of a star increases, the 
accumulate completeness converges on the ultimate completeness, Coo,j- 



3. The Probability of a Discovery by the Next LSO, V 



After i — 1 unproductive LSOs of star j, (1 — Cj_ij) is the fraction of all possible POIs 
that have not been ruled out and are still in play. The fraction of the remaining POIs that 
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Fig. 2. — The dynamic completenesses of five hypothetical LSOs of HIP 29271 at t — 10^, 
4x10^, 1.6x10''', 2.7x10''', and 3.5x10''' sec. The virgin completeness Cij = 0.46 is an accurate 
value, independent of time and computed in advance of simulations, using Eq. 5. The symbols 
C2_5,j are preceded by tildes to indicate they are approximations. Such approximations would 
be used only in the selection process of a DRM, where it may be necessary to compute 
dynamic completeness for many stars on the fiy. After the decision to observe a particular 
star, Eq. 5 would be used to produce an accurate value for the record. Note that as the 
number of LSOs of a star increases, the accumulate completeness converges on the ultimate 
completeness, Cooj- 
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the next (i* ) LSO at time t would detect is 

The probability of a discovery on the next LSO of star j is 

Vi,j = v'j}Ci,jit) , (9) 

where r]j is the Bayesian correction of the occurrence probability rj. rj'j is the probability 
that star j possesses a POI after taking into account the contrary evidence of rij previous 
unproductive LSOs with accumulated completeness Q-ij- 

Bayes's theorem states: 

where the hypothesis H is that star j has a POI; the evidence £ is the lack of a discovery so 
far; P{S\1-L) = (1 — Ci-ij) is the conditional probability of £ if H is true; PiT-L) — rj is the 
prior probability of l-L; and the marginal probability of S is 

p{£) = p{s\n)P{'H) + p{s\n)p{n) 

= {1 - Ci.ij)ri + (1 - T]) = 1 - rid.ij , (11) 

where H is the hypothesis that star j does not have a POI. [P(£|?^) = 1 and Pi^-L) — 1 — r).] 
The result for Vij is 

For the case of no prior searches {i — 1, Cqj = 0), the probability of a discovery on the first 
search is rjcij, as expected. 

We want to confirm numerically that Eq. 12 accurately estimates the probability of a 

discovery by the next LSO for random values of the parameters, for example, 77 = 0.272673, 
Cij = 0.437671, and C2,j = 0.506385, for which V2,j = 0.156789. For this verification, we 
conduct Nohs = 200, 000 independent LSOs, each involving Nq = 200, 000 possible POIs. 
Each LSO involves the following computational steps: 



1. Randomly pick the serial number of the "real" POI: n — B{ri)I{l, Nq), where B is a 
Bernoulli random deviate with probability 77, which yields the value or 1, and I is 
a uniform random deviate producing an integer in the range 1 to A^o- (If the serial 
number is zero, it means the "star" being observed does not have a POI.) 
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2. Perform a first LSO by selecting A^i = RoUND(CijAi'o) random integers in tlie range 
1 to Ai"o, where the Round function yields the closest integer to the argument. If n is 
one of these Ni integers, then go back and repeat Steps 1 and 2, because we only want 
cases where the first LSO does not make a "discovery." 

3. Perform a second LSO by randomly selecting N2 = Round(c2jA^o) integers from the 
Nq-Ni integers defined by excluding the Ni integers in Step 2 from the set of all 
integers 1 to A^o- If ^ is equal to one of these N2 integers, then we have a discovery; 
otherwise, not. 

4. Repeat Steps 1 to 3 A^obs times, and count the number of discoveries, A^disc- 

5. Compute the empirical probability, Adisc/A'obs, which was 0.156845, 0.157050, and 
0.156695 in three runs we performed using the parameters above. These values compare 
well with the theoretical value, V2,j and confirm Eq. 12. (As a benchmark, one run 
required 20,000 sec on a 3 GHz Intel Xenon processor running Mathematica 6.) 



4. Applications of to Observing Programs 

Two applications of V and Eq. f[T^ must be sharply distinguished. The first application 
is in the scheduling algorithm for real or simulated observing programs, where we use the 
discovery rate, 

= , (13) 

tcxp + OH 

as a merit function or science benefit / cost metric for optimizing the observing program for 
discovery. In this application, 77 in Eq. (12) is ?7ops. (OH is any observational overhead time 
that will be charged to the program, such as for calibration or alignment.) 

The second application of Eq. f|T2l) is in estimating the probability distribution of dis- 
coveries for simulated observing programs, as discussed below. In this application, r] = r^tmc, 
where ?7true is the "true" value. In this context, r/true is a control parameter. If r/ops 7^ ^truc, 
the number of planets discovered by an observing program may be less than optimal, because 
the scheduling algorithm may not always choose, as the next star to observe, the qualified 
star with the highest "true" value of the merit function. 

In both applications of Eq. fll2p . Cj„ij is the accumulated completeness from all LSOs 
of star j prior to the i^^, each contribution computed as accurately as desired from Eq. 
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4.1. Probability Distribution of the Number of Discoveries, pdf{m) 

No matter what star j, nor what search i of that star, the k^^ LSO in the overaU observing 
program discovers u planetary systems, where u G {0, 1} is a Bernoulli random variable with 
probability Vk, as given in Eq. f lT2|) using r/true and the indices i and j corresponding to the 
k^^ LSO. [For each star (j), multiple visits {i) are possible, so both i and j define the k^^ 
LSO.] The probability density function (pdf) of u for the k^^ LSO is 

pdUu) = (1 - Vk)S{u, 0) + VkS{u, 1) , (14) 

where = 1 for i = j, and zero otherwise, is the Kronecker delta. 

An entire observing program, consisting of ntotai LSOs, where 

stars 

JT-total = ^IT-j y (15) 

i=i 

where rij is the total number of LSOs of star j, and where nstars is the total number of 
stars observed, discovers m planets, where m G {0, 1, . . .ntotai} is the sum of ntotai Bernoulli 
random variables, each with pdfj. given by Eq. ( IT^ . Therefore, the pdf of m is the convolution 
(*) of pdfj. for all k: 

pdf{m) = pdf^{u) -k pdf^iu) -k...-k pdf^^^^Ju) , (16) 
where each successive convolution has the form 

ipdff,kJ^){n) = (1 - n)^(n) + VkJ'in - 1) . (17) 

Equation (16) offers a practical advantage for estimating the outcome of a search pro- 
gram. Starting from a single simulated observing program, it allows us to estimate the- 
oretically the pdf of the total number of discoveries. The alternative — running many full 
simulations to build up an empirical estimate of the pdf from the discovery results — is much 
less efficient. 



4.2. Estimation of 77 

The pertinent record of a real or simulated observing program is a set of ntotai data 
triplets re- indexed from i (LSO) and j (star) to k (observation): 

{cfc, Ck, Uk} = {cij, Cj_ij, Uij} , (18) 
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for 1 < A; < ritotai, where Uij = or 1 is the number of discoveries by the z*^ LSO of the j'*^ 
star. (If we assume that we stop searching after a discovery, Uij = 1 for at most one value 
of i for any j.) 

The logarithmic likelihood function C is the logarithm of the probability of the set 
{uk} = {uij} as a function of rj: 

"total "total / / \ \ 

rnukM = E lnp„, = g In (yZ^^^ + - -T^J - ""'^j ■ ^^^^ 

The maximum-likelihood estimate of the occurrence probability, E{ri), is the r^-root of the 
equation: 

^^%M = 0. (20) 
or] 

The minimum variance bound (MVB) is the inverse of the Fisher information near E{ri): 

-1 



MVB {E{ri)) 



d'C{{uk}\v) 



(21) 



We want to confirm that Eqs. (19-21) accurately estimate rj and its variance. To 
that end, we performed a numerical experiment simulating 100,000 missions of 100 LSOs, 
according to the following steps. 

1. Generated ntotai = 100 random data triplets for the left side of Eq. f|T8|) as follows: 

c = n , (22) 
C = (1 - c)7^ , (23) 

" - Ki^) • '''' 

where 7?. is a uniform random deviate on the range 0-1, B{p) is a Bernoulli random 
deviate with probability p, and rj = 0.10. 

2. Compute E{ri) and MVB(E(r/)) using Eqs. (19-21). 

For the sample of 100,000 trials, we found 

{E{r])) = 0.101 , (25) 
an = 0.039 , (26) 
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and 

(^^/MVB{E{r]))^ = 0.042 . (27) 

Equation is the mean value of r] found using Eq. fEUl) in the 100,000 simulated 
missions. Equation ( l26l) . the standard deviation of those values of r], is the empirical estimate 
of the scatter in rj determined by Eq. (120|) . The square of Eq. (126!) is the empirical variance. 
The Cramer-Rao theoretical limit on the variance of any estimator of rj is the MVB given by 
Eq. ( 12T|) . We computed the MVB for each simulated mission, and Eq. (!27|) gives the mean 
value of the square root of the MVB for the suite of 100,000 simulated missions — computed 
for direct comparison with the empirical value in Eq. (12^ . 

These results illustrate that the maximum likelihood estimator accurately recovers 1] 
from a record of the results of observing programs, and that the accuracy of this estimator 
appears to approach the Cramer- Rao limit (a^ ~ MVB). (This Monte Carlo experiment 
required 270 sec on 56 2.66 GHz Intel Xenon processors operating in parallel.) 



5. Illustrative Design Reference Missions (DRMs) 

The purpose of a design reference mission (DRM) is to gauge the science operations of 
a mission concept. To illustrate the new completeness methods introduced in this paper, 
we now describe a ministudy using simple DRMs to explore and measure of the power of 
the James Webb Space Telescope (JWST) to discover and characterize Earth- like extrasolar 



planet s using a starshade to suppress scattered starlight (ICash et al.ll2009l : ISoummer et al. 



2009al ). In this scenario, JWST and the starshade revolve in coordinated orbits around 
the second Earth-Sun Lagrange point, L2. The starshade operates on a ~70,000 km sphere 
centered on JWST. In a 3-year planet-finding mission, we assume enough propulsion to slew 
the starshade 70 times to take up new positions between JWST and target stars. We want 
the DRMs to tell us about the science, for example how many discoveries to expect if we 
optimize the observing program, assuming r^true = Vops = 0.3, say. 

Other DRM inputs include: a science strategy; a definition of POIs (same as Section [521 
with q depending on filter as given in Table 1); IWA = 0.085 arcsec, AmagQ = 26, and point- 
ing restrictions 71 = 85° (solar avoidance) and 72 = 105° (starshade bright-side avoidance); 
an input catalog of stars; exposure time calculators; typical overheads OH; and a merit 
function — in this case, the discovery rate Z — which we use to select the next star to search. 

The science strategy is simple: we perform an LSO, and if a likely POI is discovered, we 
immediately perform additional spectrophotometry to characterize the body. Such immedi- 
ate follow-up reduces the risk of a newly discovered POI becoming undetectable before it can 
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be characte rized, and avoids the difficulty of t rying to recover it with inadequate knowledge 



of its orbit (IBrown. Shaklan. &: Hunyadil 120071 ). If we discover a POI, we cease further LSOs 



of that star. If we do not find one, we move the starshade to the next target star, but we 
return to a star already searched if and when it once again offers the highest value of Z. 

The LSO is a deep image using whichever NIRCam filter offers maximum Z. We call 
this filter the preferred filter, and it varies from star to star. The possible filters are listed 
in Table 1. 

After an LSO finds a potential POI, we obtain images through the four non-preferred 
NIRCam filters for that star, and take a low-resolution spectrum with NIRSpec. We use 
exposure time calculators for NIRCam and NIRSpec, based on parameters from the instru- 
ment teams, to achieve S/N = 5 on a source of magnitude magj + Amag^ for the LSOs, and 
magj + Amag^^^^^^ for the follow-up filter photometry and spectroscopy, where mag^ is the 
apparent magnitude of star j, and Amag^^^^^^^ is the median magnitude difference between 
the star and the universe of possibly detected POIs for that star. (If Amag^^^^g^^ > AmagQ, 
we use AmagQ.) 

At this stage of the study we do not try to refine our understanding of the exposure time 
calculation beyond the current estimates by the instrument teams (Marsha Rieke et al. and 
Peter Jakobsen et al., private communication). We use standard parameters based on in- 
strument requirements. There may be better observing modes for this particular application 
(e.g., involving detector sub-arrays). 

We interpolate standard stellar magnitudes and zero points to the effective wavelengths 
of each NIRCam filter and the NIRSpec prism, starting from the VJHK magn i tudes from 



NStED (nsted.ipac.caltech.edu/) and the VJHK zero points from lLeinert et al.l ( 1l998l ). We 
used a near-infrared spectrum of Earth calculated by Sara Seager (private communication) 
to estimate the effective geometric albedo q of Earth at the wavelengths and resolving powers 
of each instrument modes. These parameters are listed in Table 1. 

As illustrated in Figure 3, each potential target star is continuously observable only for 
a limited period of time, once or twice a year, depending on its ecliptic latitude (6). The total 
time costs of observing a star — itemized in Tab le 2 — rn u st fit into a single observability period 



for that star. Some 26 of the 117 target stars iBrownl (120051 ) used in an earlier study of the 



coronagraphic Terrestrial Planet Finder ( TPF-C) are qualified according to this criterion. 
These stars constitute the input catalog for these DRMs (Table 3). 

To determine the preferred filter for LSOs, we use Eq. (5) and the procedures of Sec- 
tion |22] to calculate virgin completeness cij for each of the 26 qualified stars using samples 
of A^o = 40, 000 POIs. We do this separately for each of the five filters, because of the depen- 
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Fig. 3. — The sphere of starshade operations, centered on JWST, shown here on the vernal 
equinox. Green: permitted pointings for 71 = 85° and 72 = 105°. Red: forbidden pointings. 
As seen from above and as time passes, stars revolve on the starshade sphere around the +z 
axis (north ecliptic pole) in the clockwise direction. Depending on a star's ecliptic latitude 
b, it may be observable for one or two period s per y ear, or for the entire year {b > 85°). 
Blue: a typical TPF-C target star from|Brown| J2005I . HIP 92043, b = +43.4°). Black: other 
TPF-C stars. Cyan: the Sun, which is fixed in this L2 coordinate system (— x toward the 
Sun, +y in the direction of the Earth's orbital motion). 
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dence of q on wavelength and resolving power. Next, we use Eqs. (12-13) to compute Vij 
and Zi j, using OH = 10 hours as our estimate of the time cost of the fine alignment of the 
starshade, which we assume is incurred by each new observation with a different instrument, 
or of a different star. For each star j, we select the filter with the highest value of Zi j as 
the preferred filter for LSOs of that star (listed in Table 3). 

Next, we determine the universe of possible LSOs. For the preferred filter only, we 
continue to compute Cj>i j using Eq. (5) until the sample of Nq POIs is effectively exhausted. 
This yields 26 lists of dynamic completenesses, in sequence, one list for each of the 26 stars — 
some 2,075 values of Cij in all. Again using Eqs. (12-13), we convert these hsts of Cij into 
a full list of possible LSOs in the form of vectors {HIPfc, ik, Vk, 2k}, where 1 < k < 2075 
is the index for LSOs introduced in Section 14. 1^ and the items are the Hipparcos number 
of the star, the number of this visit to that star, and the discovery probability and rate for 
that visit. 

At the start of a DRM, the prioritized observing program is the list of 2,075 LSO vectors 
sorted in descending order of Zk- (Table 4 lists the top 80 LSOs for the current illustration.) 
We expect to perform only 70 LSOs — but we do not know which ones. How far down the 
list a DRM reaches is determined by the random discoveries as the DRM unfolds. That is, 
we determine the outcome of each LSO in turn — discovery, yes or no, with possible follow- 
on observations and alternative time costs in Table 2 — by interrogating a Bernoulli random 
deviate with probability V, and interpreting "1" as a discovery and "0" as no discovery. 
Each discovery deletes from the observing program all the pending LSOs of a star — meaning 
those with i greater than the visit number i of the LSO that produced the discovery. These 
deletions promote the lower priority LSOs of other stars into higher positions on the list. 
To investigate this behavior and its ramifications, we conduct a Monte Carlo experiment of 
500,000 DRMs. 

Table 4 provides an example DRM in the form of the LSOs actually executed in one 
DRM run, and the discoveries actually made — five in this case. We use it to illustrate the 
method of deriving the theoretical probability distribution of m from the probabilities Vk of 
a single DRM. We collect the 70 values of Vk^^^^ for 1 < A;drm < 70 in Table 4 and follow the 
recipe in Section UT] to obtain the theoretical pdf{m) represented by the dots in Figure 4. 
For comparison, the histogram in Figure 4 shows the empirical pdf{m) based on the actual 
values of m from 500,000 DRMs. 

Table 5 compares the means, standard deviations, and standard deviations of the 
means — (m), am, and cr^m) — for the results ensuing from Table 4, which encompass both 
the example DRM and the 500,000 DRMs computed from the same suite of potential or 
actual LSOs (1 < A; < 2075 or 1 < /cdrm < 70). The highest precision is achieved by the 
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Fig. 4. — Comparison of the theoretical and empirical pdf{'m). Dots: the theoretical pdf 
estimated from the 70 values of Vk of the actual LSOs in the single example DRM in Table 4. 
Histogram: the empirical pdf from the actual results for m in 500,000 DRMs. (The histogram 
has been normahzed to 1.0.) 
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mean of the 500,000 individual theoretical results for (m). We see that our new method 
produces accurate theoretical estimates of (m) and am from a single DRM and avoids the 
onerous alternative, which is to conduct a large number of DRMs to obtain empirical results. 

The 500,000 DRMs ensuing from the prioritized LSOs in Table 4 offer additional insights 
into the science operations of the JWST starshade mission. The most likely number of unique 
stars searched is 23 (ranges 22-25). The median number of visits per star is 3. The total 
JWST observing time used by a DRM is 7.6 dz 1.3 x 10^ sec, and the total number of LSOs is 
70 (here, this is the limiting parameter). Note that the total observing time corresponds to 
only ~7% of the total JWST observing time. The most likely probability of zero discoveries 
is 0.004 (ranges 0.003 to 0.006). 

Table 6 gives the data triplets defined in Eq. f ITS]) . which we use to estimate r^true = 
0.30 ± 0.12, where the quoted error is the square root of the MVB. 

A last note on this ministudy. We have not treated the recovery times of discussed 
in Section 12.31 because these example DRMs are "dilute," meaning the LSOs account for 
only a small fraction of all the observations of JWST. In a "dense" DRM, with every LSO 
competing to be the next observation of the telescope, recovery time is important. Here, 
however, we can assume that a buffer of at least ~10^ sec can conveniently pad the time 
between any two visits of the same star, which ensures adequate Cij recovery. 

6. Summary 

In this paper, we have extended completeness-based metrics and algorithms for direct 
exoplanet searches to include multiple visits, estimating the probability distribution of search 
results, and estimating the occurrence rate of extrasolar planets. These extensions open the 
way for improved scheduling decisions, more realistic expectations for candidate instruments 
and missions, and enhanced science returns. 

Preliminary DRMs for a starshade mission to enable Earth-like planet searches with 
JWST suggest a viable program for < 10^ sec of observing time. About five discoveries and 
spectral characterizations are expected if 77 = 0.3. Soon, we hope, the Kepler mission will 
provide a first estimate of the true value of r]. 

The authors thank Jay Anderson, Massimo Robberto, Elizabeth Barker, Marsha Rieke, 
and Doug Kelly for help with the NIRCam parameters; Jeff Valenti, Jason Tumlinson, Peter 
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Table 1. Spectrophotometric parameters for calculating completenesses and exposure 

times. 







Nominal A 


Resolving 


Zero point 


Earth geometric 


Instrument 


Mode 


(nm) 


power 


(Jy) 


albedo (q) 




F070W 


700 


4 


3043 


0.232 




F115W 


1150 


4 


1766 


0.187 


NIRCam 


F140M 


1400 


10 


1324 


0.021 




F150W 


1500 


4 


1188 


0.103 




F162M 


1625 


10 


1045 


0.179 


NIRSpec 


])rism 


1150 


31.7 


1766 


0.260 



Note. — Zero points interpolated from values for VJHK filters (Leinert 1997). 
We adjusted the Earth's effective albedo for prism spectroscopy at 1150 nm. 
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Table 2. Observing sequence and typical time costs for a single visit of the starshade to a 

target star. 



Step 


Activity 


Observing time cost 


Clock time cost 


1 


Final alignment of JWST, starshade, and 


10 hours 


10 hours 




target star 






2 


LSO through preferred NIRCam filter 


4.8 X 10^ sec 


4.8 X 10^ sec 


3 


Analyze data for discovery 




7 days 


4 


Final ahgnment of JWST, starshade, and 


10 hours 


10 hours 




target star 






5 


Images through four non-preferred NIR- 


1.4 X 10^ sec 


1.4 X 10^ sec 




Cam filters 






6 


Analyze data to prepare for spectroscopy 




7 days 


7 


Final alignment of JWST, starshade, and 


10 hours 


10 hours 




target star 






8 


Spectrum using NIRSpec 


1.9 X 10^ sec 


1.9 X 10^ sec 




Total with discovery 


4.9 X 10^ sec 


1.7 X 10*' sec 




Total without discovery 


8.4 X 10^ sec 


6.9 X 10^ sec 



Note. — Overhead costs (1, 3, 4, 6, 7) are fixed. Exposure times (2, 5, 8) are median 
values for the 26 stars in the input catalog. With no discovery in step 3, only steps 1-3 are 
executed during a visit to a target star. With a discovery, all 8 steps are executed. During 
the analysis steps (3, 6), JWST conducts observations for other programs; only the starshade 
remains aligned with the target star for this program. 



Table 3. The input catalog of target stars, ranked in descending order of the discovery rate on the first LSO {2i). 



Max 



Disc. 








Distance 






T Qf^ 


A /Trt .... 

Max 


cont. 










rate 




lype 


L 




o( ) 


nlter 


''exp 


time 


obs. 


Cl 




ri 




1 


71681 


Kl V 


0.61 


1.35 


-43 


F070W 


3.2 


6.12 


6.38 


0.8 


1 


0.24 


-5.19 


2 


8102 


G8 V 


0.47 


3.65 


-25 


F115W 


4.03 


6.13 


6.29 


0.76 


1 


0.23 


-5.31 


3 


71683 


G2 V 


2.2 


1.35 


-43 


F070W 


2.43 


6.12 


6.38 


0.55 


1 


0.17 


-5.34 


4 


3821 


GO V 


1.2 


5.95 


47 


F115W 


3.94 


6.14 


6.41 


0.59 


1 


0.18 


-5.4 


5 


99240 


G6/8 IV 


1.5 


6.11 


-45 


F115W 


4.12 


6.16 


6.4 


0.56 


1 


0.17 


-5.47 


6 


108870 


K4/5 V 


0.2 


3.63 


-41 


F150W 


4.48 


6.13 


6.37 


0.65 


1 


0.19 


-5.53 


7 


22449 


F6 V 


2.6 


8.03 


-15 


F115W 


3.88 


6.15 


6.26 


0.39 


0.96 


0.12 


-5.57 


8 


19849 


KO/1 V 


0.41 


5.04 


-28 


F115W 


4.69 


6.16 


6.3 


0.64 


1 


0.19 


-5.64 


9 


15510 


G8 III 


0.71 


6.06 


-58 


F115W 


4.68 


6.18 


6.53 


0.63 


1 


0.19 


-5.65 


10 


2021 


Gl IV 


3.9 


7.47 


-65 


F070W 


4.2 


6.15 


6.63 


0.34 


0.92 


0.1 


-5.7 


11 


27072 


F6.5 V 


2.3 


8.97 


-46 


F115W 


4.47 


6.25 


6.4 


0.41 


0.98 


0.12 


-5.73 


12 


1599 


GO V 


1.2 


8.59 


-58 


F115W 


4.71 


6.25 


6.53 


0.5 


1 


0.15 


-5.76 


13 


64394 


GO 


1.3 


9.15 


33 


F115W 


4.81 


6.3 


6.32 


0.47 


1 


0.14 


-5.85 


14 


57757 


F8 


3.4 


10.9 


0.69 


F115W 


4.33 


6.23 


6.24 


0.27 


0.88 


0.08 


-5.85 


15 


12777 


F8 


2.2 


11.2 


32 


F115W 


4.64 


6.31 


6.32 


0.36 


0.99 


0.11 


-5.87 


16 


96100 


G9 V 


0.41 


5.77 


81 


F115W 


4.99 


6.2 


7.42 


0.56 


0.97 


0.17 


-5.9 


17 


105858 


F7 V 


1.4 


9.22 


-47 


F115W 


4.91 


6.33 


6.41 


0.47 


1 


0.14 


-5.92 


18 


73184 


K4 V 


0.26 


5.91 


-4.4 


F150W 


5.05 


6.2 


6.25 


0.34 


0.76 


0.1 


-6.17 


19 


70497 


F8 


3.9 


14.6 


60 


F115W 


4.75 


6.41 


6.56 


0.17 


0.81 


0.051 


-6.26 


20 


23693 


F6/7 V 


1.4 


11.7 


-79 


F115W 


5.18 


6.53 


7.42 


0.33 


0.91 


0.098 


-6.29 


21 


77952 


FO III/IV 


9.6 


12.3 


-42 


F070W 


4.3 


6.23 


6.38 


0.073 


0.4 


0.022 


-6.41 



Table 3 — Continued 



Max 



Disc. 








Distance 




LSO 


LSO 


Max 


cont. 










rate 


HIP 


Type 


L 


(pc) 


K) 


filter 


^exp 


time 


obs. 


Cl 


C'oo 






22 


29271 


G6 V 


0.83 


10.2 


-82 


F115W 


5.41 


6.53 


7.43 


0.32 


0.81 


0.095 


-6.49 


23 


114622 


K3 V 


0.28 


6.53 


55 


F150W 


5.34 


6.27 


6.49 


0.27 


0.68 


0.081 


-6.5 


24 


50954 


F2/3 IV/V 


5 


16.2 


-68 


F115W 


4.98 


6.57 


6.7 


0.11 


0.68 


0.032 


-6.61 


25 


40702 


F5 V 


6.8 


19.4 


-75 


F070W 


5.29 


6.66 


7.41 


0.076 


0.61 


0.023 


-7.01 


26 


86611 


F5 


5.5 


22 


81 


F070AA' 


5.6 


6.93 


7.16 


0.018 


0.67 


0.011 


-7.48 



Note. — Times in log seconds. "Max time" is the time cost of the full observing sequence in Table 1. "Max 
cont. obs." is the maximum continuous observing time for ecliptic latitude h for 71 = 85° and 72 = 105°. Vi is 
the probability of a discovery on the first LSO, from Eq. (12), and Zi is Vi divided by the sum of LSO texp and 
10 hours for alignment, from Eq. (13). The value of Zi is given in log discoveries per sec. 
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Table 4. Record of a typical DRM. 



k 




HIP 


i 


Vh 

' K 




Discovery 


k 


"'Uri.lvi 


HIP 


i 


' k 




Discovery 


1 


1 


71681 


1 


0.24 


-5.19 


no 


41 


37 


57757 


3 


0.039 


-6.17 


no 


2 


2 


8102 


1 


0.23 


-5.31 


yes 


42 


38 


105858 


2 


0.075 


-6.19 


no 


3 


3 


71683 


1 


0.17 


-5.34 


no 


43 




71683 


4 


0.022 


-6.22 




4 


4 


3821 


1 


0.18 


-5.40 


no 


44 


39 


12777 


3 


0.046 


-6.24 


no 


5 


5 


99240 


1 


0.17 


-5.47 


yes 


45 


40 


2021 


4 


0.029 


-6.26 


no 


6 


6 


108870 


1 


0.19 


-5.53 


no 


46 


41 


70497 


1 


0.051 


-6.26 


no 


7 


7 


22449 


1 


0.12 


-5.57 


no 


47 


42 


23693 


1 


0.098 


-6.29 


no 


8 


8 


19849 


1 


0.19 


-5.64 


no 


48 


43 


57757 


4 


0.029 


-6.30 


no 


9 


9 


15510 


1 


0.19 


-5.65 


no 


49 


44 


1599 


3 


0.043 


-6.31 


no 


10 


10 


71683 


2 


0.079 


-5.66 


yes 


50 


45 


96100 


2 


0.065 


-6.31 


no 


11 


11 


2021 


1 


0.10 


-5.70 


no 


51 


46 


108870 


3 


0.032 


-6.32 


no 


12 


12 


27072 


1 


0.12 


-5.73 


no 


52 




8102 


3 


0.022 


-6.32 




13 


13 


1599 


1 


0.15 


-5.76 


no 


53 


47 


27072 


4 


0.030 


-6.34 


no 


14 


14 


3821 


2 


0.077 


-5.77 


no 


54 


48 


4394 


3 


0.045 


-6.35 


no 


15 


15 


22449 


2 


0.070 


-5.79 


no 


55 


49 


71681 


3 


0.017 


-6.35 


no 


16 




99240 


2 


0.079 


-5.79 




56 




99240 


4 


0.022 


-6.35 




17 


16 


71681 


2 


0.054 


-5.84 


no 


57 


50 


22449 


5 


0.019 


-6.36 


no 


18 


17 


64394 


1 


0.14 


-5.85 


no 


58 


51 


70497 


2 


0.040 


-6.37 


no 


19 


18 


57757 


1 


0.080 


-5.85 


no 


59 




3821 


4 


0.019 


-6.37 




20 


19 


12777 


1 


0.11 


-5.87 


no 


60 


52 


15510 


3 


0.035 


-6.38 


no 


21 




8102 


2 


0.059 


-5.90 




61 


53 


12777 


4 


0.033 


-6.39 


no 


22 


20 


96100 


1 


0.17 


-5.90 


no 


62 


54 


2021 


5 


0.020 


-6.41 


yes 


23 


21 


2021 


2 


0.063 


-5.91 


no 


63 


55 


77952 


1 


0.022 


-6.41 


no 


24 


22 


105858 


1 


0.14 


-5.92 


no 


64 


56 


105858 


3 


0.045 


-6.42 


no 


25 


23 


27072 


2 


0.073 


-5.96 


no 


65 


57 


19849 


3 


0.032 


-6.42 


no 


26 




71683 


3 


0.040 


-5.96 




66 


58 


57757 


5 


0.020 


-6.45 


no 


27 


24 


22449 


3 


0.044 


-6.00 


no 


67 


59 


70497 


3 


0.031 


-6.48 


no 


28 


25 


108870 


2 


0.065 


-6.01 


no 


68 


60 


29271 


1 


0.095 


-6.49 


no 


29 


26 


57757 


2 


0.055 


-6.02 


no 


69 


61 


73184 


2 


0.048 


-6.49 


no 


30 


27 


1599 


2 


0.075 


-6.06 


no 


70 


62 


22449 


6 


0.014 


-6.49 


no 


31 


28 


12777 


2 


0.069 


-6.07 


no 


71 




71683 


5 


0.012 


-6.50 




32 


29 


15510 


2 


0.071 


-6.07 


no 


72 


63 


23693 


2 


0.060 


-6.50 


no 


33 


30 


3821 


3 


0.038 


-6.08 


yes 


73 


64 


114622 


1 


0.081 


-6.50 


no 


34 




99240 


3 


0.041 


-6.08 




74 


65 


27072 


5 


0.021 


-6.50 


no 


35 


31 


2021 


3 


0.042 


-6.09 


no 


75 


66 


1599 


4 


0.026 


-6.52 


no 


36 


32 


19849 


2 


0.068 


-6.09 


no 


76 


67 


77952 


2 


0.017 


-6.52 


no 


37 


33 


64394 


2 


0.074 


-6.13 


no 


77 


68 


12777 


5 


0.022 


-6.56 


no 


38 


34 


27072 


3 


0.047 


-6.15 


no 


78 


69 


57757 


6 


0.016 


-6.56 


no 


39 


35 


22449 


4 


0.030 


-6.16 


no 


79 


70 


64394 


4 


0.028 


-6.56 


no 


40 


36 


73184 


1 


0.10 


-6.17 


no 


80 




2021 


6 


0.014 


-6.57 





Note. — k: order of top 80 LSOs in terms of discovery rate Z; feoRM: order of 70 LSOs executed in this example DRM; 
i: number of this LSO of this star; Vk- probability of discovering a POI in this LSO; Z}^: log discovery rate of this LSO; 
ellipses indicate that this LSO was not peformed in this DRM due to a previous discovery. 
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Table 5. DRM results for (m), am, and (T(m,) for the expected number of POIs discovered. 





Empirical 


Theoretical 




(m) am CT{m) 


(m) am cr(m) 


Example DRM 


5 


5.21 2.15 


500,000 DRMs 


5.132 2.06 0.003 


5.129 0.09 0.0001 



Note. — Empirical results are based on counting the number 
of discoveries in Monte Carlo experiments where discoveries were 
decided by Bernoulli random deviates. Theoretical results use 
the methods of Section UT] to estimate pdf{m) from the individ- 
ual LSO probabihties (Vk), and the results are — or are derived 
from — that distribution. (The highest precision is achieved by 
the mean of the 500,000 individual theoretical results for (m).) 
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Table 6. Data triplets for estimation of 77 from results of the typical DRM. 



k 




Cfe 




k 








k 




Cfe 


Uk 


1 


0.8038 








25 


0.1759 


0.6484 





48 


0.1203 


0.6833 





2 


0.7574 





1 


26 


0.1681 


0.2672 





49 


0.04000 


0.9398 





3 


0.5537 








27 


0.2123 


0.5008 





50 


0.04858 


0.7929 





4 


0.5950 








28 


0.2043 


0.3599 





51 


0.1264 


0.1697 





5 


0.5560 





1 


29 


0.1932 


0.6266 





52 


0.08750 


0.8198 





6 


0.6484 








30 


0.09530 


0.8049 


1 


53 


0.08603 


0.6926 





7 


0.3886 








31 


0.1171 


0.5334 





54 


0.05295 


0.7275 


1 


8 


0.6446 








32 


0.1827 


0.6446 





55 


0.07265 








9 


0.6266 








33 


0.2121 


0.4711 





56 


0.1190 


0.6811 





10 


0.2199 


0.5537 


1 


34 


0.1264 


0.6184 





57 


0.07973 


0.8273 





11 


0.3442 








35 


0.07820 


0.7147 





58 


0.05505 


0.6267 





12 


0.4061 








36 


0.3361 








59 


09337 


2960 





13 


0.5008 








37 


0.1120 


0.4353 





60 


0.3181 








14 


0.2099 


0.5950 





38 


0.2154 


0.4657 





61 


0.1448 


0.3361 





15 


0.2067 


0.3886 





39 


0.1284 


0.5642 





62 


0.03500 


0.8414 





16 


0.1360 


0.8038 





40 


0.07705 


0.6505 





63 


0.1806 


0.3261 





17 


0.4711 








41 


0.1697 








64 


0.2701 








18 


0.2672 








42 


0.3261 








65 


0.05180 


0.8227 





19 


0.3599 








43 


0.07945 


0.5473 





66 


0.06625 


0.8257 





20 


0.5590 








44 


0.1126 


0.7131 





67 


0.05485 


0.07265 





21 


0.1893 


0.3442 





45 


0.1813 


0.5590 





68 


0.05650 


0.7786 





22 


0.4657 








46 


0.07942 


0.8242 





69 


0.04180 


0.6818 





23 


0.2123 


0.4061 





47 


0.07785 


0.7449 





70 


0.06988 


0.8035 





24 


0.1193 


0.5953 






















Note. — k here is ^drm in Table 4. Vk in Table 4 can be recovered from Ck and Ck here 
using Eqs. (12 h 18) and t) = r/ops = 0.3. 



